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Abstract 

Accumulating evidence suggests that many tumors have a hierarchical organization, with the bulk of the tumor composed 
of relatively differentiated short-lived progenitor cells that are maintained by a small population of undifferentiated long- 
lived cancer stem cells. It is unclear, however, whether cancer stem cells originate from normal stem cells or from 
dedifferentiated progenitor cells. To address this, we mathematically modeled the effect of dedifferentiation on 
carcinogenesis. We considered a hybrid stochastic-deterministic model of mutation accumulation in both stem cells and 
progenitors, including dedifferentiation of progenitor cells to a stem cell-like state. We performed exact computer 
simulations of the emergence of tumor subpopulations with two mutations, and we derived semi-analytical estimates for 
the waiting time distribution to fixation. Our results suggest that dedifferentiation may play an important role in 
carcinogenesis, depending on how stem cell homeostasis is maintained. If the stem cell population size is held strictly 
constant (due to all divisions being asymmetric), we found that dedifferentiation acts like a positive selective force in the 
stem cell population and thus speeds carcinogenesis. If the stem cell population size is allowed to vary stochastically with 
density-dependent reproduction rates (allowing both symmetric and asymmetric divisions), we found that dedifferentiation 
beyond a critical threshold leads to exponential growth of the stem cell population. Thus, dedifferentiation may play a 
crucial role, the common modeling assumption of constant stem cell population size may not be adequate, and further 
progress in understanding carcinogenesis demands a more detailed mechanistic understanding of stem cell homeostasis. 
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Introduction 

Most tissues consist of three classes of cells: stem cells, transit- 
amplifying progenitor cells, and differentiated cells. Multicellular 
organisms require a tight control of cell division to ensure a proper 
balance between these different cell populations. The cancer stem 
cell (CSC) hypothesis states that tumors are also hierarchically 
organized, with a small sub-population of cancer cells driving 
cancer growth [1]. Individual cell tracing studies of tumor 
development strongly support the cancer stem cell hypothesis in 
many (but not all) types of cancer [2,3], and identifying these cells 
in tissues is an ongoing goal in cancer research. Lineage studies 
find that malignant tumors contain more cancer stem cells 
compared to benign tumors and that cancers gradually lose their 
tissue-like hierarchical organization as they evolve from benign to 
malignant [2]. 

Cells escape proliferation control after acquiring a series of 
mutations in a multi-step process [4], While some cancers may 
require only a few mutations [5], the number of required (driver) 
mutations in solid cancers is larger, with up to twenty driver 
mutations being required [6] . In order to accumulate this critical 
number of mutations during a lifetime, cells either have to be long- 
lived or the mutation rate has to be large [7] . Stem cells have been 
proposed to be likely candidates for the initial cell of mutation due 
to their long lifetime and sustained self-renewal capacity [1]. In 



addition to their long life span, stem cells are able to generate full 
lineages of differentiated cells, thereby perpetuating mutations 
through clonal expansion. Given known division and mutation 
rates, theoretical studies have argued that the necessary number of 
mutations for carcinogenesis cannot be obtained in the stem cell 
population on a reasonable time scale without assuming either 
significant selective advantage or elevated mutation rates [4,7]. 
However, there is conflicting evidence as to how early in tumor 
development cancers acquire an elevated mutation rate [8,9] and 
several cancer genome sequencing studies have estimated muta- 
tion rates during cancer initiation to be normal for some types of 
cancer [10-12]. 

Although a stem cell may sustain the first oncogenic hit, 
subsequent alterations required for development of CSCs can 
occur in descendent progenitor cells [13]. Dysregulation of 
pathways involved in stem cell self-renewal may lead to progenitor 
cells acquiring a stem cell-like phenotype. It remains an open 
question whether cancer stem cells originate from stem cells that 
escape homeostasis or from dedifferentiated progenitor cells that 
acquire infinite proliferating potential [14]. There is significant 
evidence that dedifferentiation can play a role in establishment of 
certain cancers [15—17]. For example, cell sorting has demon- 
strated that stem-like cells can arise de novo from non-stem-like 
cancer cells in in vitro breast cancer cell lines [18,19]. In the 
hematopoietic system, it has been shown that leukemic stem cells 
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Author Summary 

Recent evidence suggests that, like many normal tissues, 
many cancers are maintained by a small population of 
immortal stem cells that divide indefinitely to produce 
many differentiated cells. Cancer stem cells may come 
directly from mutation of normal stem cells, but this route 
demands high mutation rates, because there are few 
normal stem cells. There are, however, many differentiated 
cells, and mutations can cause such cells to "dedifferen- 
tiate" into a stem-like state. We used mathematical 
modeling to study the effects of dedifferentiation on the 
time to cancer onset. We found that the effect of 
dedifferentiation depends critically on how stem cell 
numbers are controlled by the body. If homeostasis is 
very tight (due to all divisions being asymmetric), then 
dedifferentiation has little effect, but if homeostatic 
control is looser (allowing both symmetric and asymmetric 
divisions), then dedifferentiation can dramatically hasten 
cancer onset and lead to exponential growth of the cancer 
stem cell population. Our results suggest that dedifferen- 
tiation may be a very important factor in cancer and that 
more study of dedifferentiation and stem cell control is 
necessary to understand and prevent cancer onset. 

can be generated from committed progenitor cells that acquire 
stem cell-like behavior [20]. It has been suggested that acute 
myeloid leukemia (AML) is a progenitor disease, where a 
progenitor acquires abnormal self-renewal potential and "dedif- 
ferentiates" to a stem cell-like state [21,22]. Other myeloid 
leukemias such as CML (chronic myeloid leukemia) are thought of 
as stem-cell diseases [23]. However, although a hematopoietic 
stem cell is thought to be the cell of origin in the early phases of 
CML, in patients with CML blast crisis, granulocyte-macrophage 
progenitors are thought to acquire self-renewal capacity through a 
/?-catenin mutation and emerge as the probable CSCs [24]. Using 
mathematical modeling to investigate the likelihood of mutation 
occurring in a progenitor versus a stem cell is a continuing line of 
investigation [25] . We treat the probability of a mutant progenitor 
cell acquiring stem cell-like state as a "dedifferentiation" rate, and 
we study how this parameter influences the time to carcinogenesis. 
We are primarily interested in whether dedifferentiation can speed 
up the time to tumor development in hierarchically organized 
cancers and in what rates of dedifferentiation are necessary for a 
noticeable effect. 

Prior Related Mathematical Modeling 

Certain aspects of the cancer stem cell hypothesis have 
previously been addressed by mathematical models. It has been 
shown that having a hierarchical tissue design, where a small 
population of stem cells maintains a transient population of 
differentiating cells, may slow the accumulation of mutations and 
protect against cancer [26-28]. The question of whether genetic 
instability (resulting in hyperactive mutation rate) is an early or 
later event in mutation acquisition leading to cancer has been 
addressed by several groups (see [4] for review). Most mathemat- 
ical models find that the onset of genetic instability should be an 
early event, if at least some of the mutations are neutral. However, 
sequencing suggests that the mutator phenotype is expressed 
relatively late in cancer progression [9]. 

Stem cell populations are typically small. Hence, the dynamics 
of mutant cells in the stem cell population are highly sensitive to 
stochastic fluctuations. A tumor begins with a single mutated cell, 
so there is a substantial chance of mutant extinction due to 



random events. Genetic drift and stochastic clonal extinction in 
stem cell lineages have been experimentally demonstrated for both 
normal tissue stem cells [29-31] and cancer stem cells [2] in 
several tissue types. Consequently, a deterministic model of 
mutation acquisition in stem cells will significantly underestimate 
the time to cancer establishment [32]. Many models of mutation 
acquisition use a stochastic approach and are concerned with 
calculating time to emergence or fixation (or when the number of 
mutant cells reaches some threshold value used in diagnosis) of a 
mutant cell with fitness r = 1 + s in a population of size N sc . 

The waiting time for cancer is often defined as the time until a 
particular number of mutation events have occurred in at least one 
cell. Iwasa et al. [33] considered a two-stage Moran model and 
described conditions under which "stochastic tunneling" can 
occur. (In this phenomenon, cells with two mutations reach 
fixation before cells with one mutation reach fixation.) Durrett et 
al. [34] obtained asymptotic estimates of waiting times until a cell 
with i mutations first appears under the assumption of neutrality 
(s = 0). These models typically consider a fixed population size 
[5,23,25,35-39]. The fixed population assumption is supposed to 
reflect homeostasis in the stem cell population, though how 
homeostasis is achieved is typically not addressed. Although the 
Moran model captures the stochastic nature of mutation 
acquisition, this type of model is not capable of describing 
mutations that change the stem cell division pattern and result in 
possible expansion of the stem cell pool, which in turn leads to 
tumor growth. Some recent models also consider mutation 
accumulation in exponentially growing cell populations [40-43]. 
Beerenwinkel et al. [6] used the Wright-Fisher model with 
exponentially growing population size to look at the effect of 
selection on the waiting time to cancer, and they predicted that the 
observed genetic diversity of colorectal cancer genomes can arise 
under a normal mutation rate (taken to be «=10~ 7 per cell 
division) if the average selective advantage per mutation is on the 
order of 1%. Similar calculations using a discrete branching 
process found 5 = 0.4% given w=10~ 5 [40]. Note that increased 
mutation rates due to genetic instability would allow even smaller 
selective advantages during tumorigenesis, but neutral mutants 
ls = 0) result in waiting times that are too long compared with 
disease incidence. Other groups have also concluded that for 
normal mutation rates and neutral mutants, mutations in multiple 
genes in acquired hematopoietic disorders are most likely very rare 
events, as acquisition of multiple mutations typically requires 
development times that are too long compared to disease 
incidence [36]. 

Spencer et al. [44] and Ashkenazi et al. [45] have focused on the 
sequential order of mutations associated with increased rate of 
proliferation, decreased rate of death, increased mutation rate, and 
other hallmarks of cancer that must accumulate before emergence 
of cancer. The sequence of mutations with the shortest waiting 
time to getting all the necessary mutations is considered the most 
likely mutational pathway [25,44]. However, these models do not 
consider the possibility that dedifferentiation of progenitor cells 
can affect the time to multiple mutation acquisition. 

The dividing progenitor cell population has previously been 
described by multi-compartment ODE models, with cells moving 
between compartments as they age [45-47]. Note that in these 
models the exact number of different stages of differentiation is 
ambiguous and does not exactly correspond to mitotic events, as 
cells may undergo more than one division in each compartment 
stage [46] . Most of these models of age-structured cell populations 
assume a stem cell proliferation rate that is dependent on the total 
number of cells and thus incorporate negative feedback as a means 
of achieving homeostasis [48,49] . These deterministic models have 
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focused on mechanisms that could regulate cell numbers that are 
necessary for homeostasis and efficient repopulation. We use a 
similar mathematical approach to model the progenitor popula- 
tion as [49] , but we couple it to stochastic dynamics in the stem cell 
compartment. 

Upon division a stem cell can produce zero, one, or two stem 
cells with probabilities op, oy, and as, respectively (Fig. 1A). The 
mean number of stem cell offspring is given by a. a -\- las. If 
symmetric divisions are permitted, the stem cell population can be 
described by a branching process with the expected number of 
cells at time t given by (a a +2as)'. However, a branching process 
either goes extinct or undergoes exponential growth, and thus it 
cannot capture stem cell dynamics at equilibrium. One solution is 
to use a conditional branching process [50], where the probabil- 
ities for a branching process are conditioned to the total 
population size remaining constant by an unspecified sampling 
mechanism (i.e., assuming that the stem cell population remains in 
homeostasis). Some theoretical studies have previously considered 
the impact of the asymmetry of cell division on stem cell dynamics. 
However, these stochastic models all assumed a fixed stem cell 
population size, either through a variant of the Moran process 
[35,51] or conditional branching process [39]. We utihze a 
different approach to get a time-varying but bounded stem cell 
population size in our models. 

Our Modeling Approach 

We use mathematical modeling to study how the possibility of 
"dedifferentiation" of mutant progenitor cells into a stem cell-like 
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Figure 1. Schematic representation of our model. (A) Stochastic 
model for stem cell division. A stem cell can produce zero, one, or two 
stem cells with probabilities c*o, ola, and as, respectively. The mean 
number of offspring is given by m = a. A +2oes. (B) Stem cells serve as an 
input to the proliferating progenitor population, and the progenitor 
population feeds back to the stem cell pool via dedifferentiation. (C) 
Mutation occurs with rate u during division and can affect cells both in 
the stem cell and progenitor pools. Blue circles represent wild-type 
cells, red circles cells with one mutation, and green circles cells with two 
mutations. (D) The sequence of models explored in this paper. 
doi:1 0.1 371 /journal.pcbi.1 003481. g001 



state affects the waiting time to carcinogenesis. Dividing progen- 
itor cells have large growing populations, so we use a deterministic 
model to describe their evolutionary dynamics. For stem cell 
populations, stochastic effects are important, because the prolif- 
erating stem cell population is typically small. We use a stochastic 
model for stem cell dynamics as a boundary condition to the PDE 
governing differentiated cell expansion (Fig. IB and C.) There is 
also feedback from the deterministic progenitor population to the 
stochastic stem cell population as a rate of "dedifferentiation". 

To assess the effect of dedifferentiation on time to carcinogen- 
esis, we consider models for stem cell dynamics with both fixed 
and variable stem cell numbers (Fig. ID). The main questions we 
address are: 

1 . What is estimated time to carcinogenesis (acquisition of M 
mutations) in stem cell-driven cancers if dedifferentiation from 
the progenitor population is allowed? 

2. What magnitude of dedifferentiation rate is needed to 
significantly shorten the time to cancer acquisition? Will 
dedifferentiation still change the waiting time to cancer if 
homeostasis in the stem cell population is maintained 
(population size remains constant) or does homeostasis need 
to be lost? 

3. What is the effect of symmetric division of stem cells, which 
leads to a non-constant stem cell population size? Do stochastic 
fluctuations in the size of the stem cell pool shorten the time to 
malignancy compared to a constant stem cell population size? 

Our general compartment model can be applied to different 
tissues, such as colonic crypts, mammary cells, and hematopoiesis. 

Models 

We make the following assumptions: 

1 . The population of progenitor cells is large enough that 
maturation can be treated as a continuous variable. Discretiz- 
ing the progenitor cell population based on the number of 
divisions a cell has completed, we obtain an age-structured 
partial-differential-equation model for the number of differen- 
tiated cells of age a at time /. We assume that when progenitor 
transit-amplifying cells carrying i mutations divide, they 
produce progenitor cells of the same maturity stage, obtaining 
the linear PDE 



dp -, da dpi 
dt dt 8a 



--{a(a)-^(a))pi. 



(1) 



where Pi(a,i) is the progenitor cell density at age a and time t, 
a(d) is the age-dependent proliferation rate, and /((a) is the 
age-dependent mortality rate. We assume that the rate of 
da 

maturation — does not depend on age a and, without loss of 

generality, set it equal to 1. Similar age-structured population 
equations have been previously studied, with focus on the 
regulatory feedback mechanisms that are necessary for homeo- 
stasis and structural stability of the steady state solution [47,49]. 
2. We assume that there is a separate stem cell population that 
gives rise to newly born differentiated cells Pi(0,t) that serves as 
a boundary conditions to the PDE system in Eq. (1). We treat 
this population stochastically. We consider two different models 
for the stem cell population: one where the total stem cell 
population is fixed and a variant that allows a time-varying but 
bounded stem cell population size (Fig. ID). 
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3. We assume neutral fitness of mutant stem cells, with the 
proliferation advantage of the mutant phenotype appearing 
only in the progenitor stage, in line with what is known for 
some cancers of the hematopoietic system [52,53]. 

4. We require M = 2 mutations to appear in the progenitor 
population before dedifferentiation to a stem-cell like state is 
possible, because sequencing of acute myeloid leukemia 
genomes suggest that there are two driver mutations present 
[54]. Additional justification for requiring M = 2 mutations is 
considered in the Discussion. 

5. Due to lack of data on dedifferentiation capacities, we assume 
progenitor cells of all maturity stages have an equal probability 
of dedifferentiating. 

Progenitor Cells 

Extending Eq. (1) to account for mutations between multiple 
subpopulations of progenitor cells (Fig. 1C) we obtain 

^^((l-^ota)-^, (2a) 



-jr~ + -p- =((^-u*)(Ti(a)-n l (a))pi+u*(To(a)po, ■ ■ • (2b) 



To mimic a maturity switch for cellular proliferation and death, 
we took the proliferation and death rates of differentiated cells per 
unit time to be 

f(«)= ^(l-tanh(p A («-coj))), and (5a) 



Ka)= -(l+tanhO/a-cOd))). (5b) 

Here b and d are the maximal proliferation and death/removal 
rates of progenitor cells. The age at which the proliferation switch 
occurs (i.e., half the progenitor cells stop dividing) is given by cu/,, 
and the steepness of the proliferation switch is determined by p h . 
Similarly, the age at which half the cells begin to undergo 
apoptosis is given by Wb, and the steepness of the death switch is 
controlled by p d . If ojf, < coj, then differentiated cells between the 
ages of co/, and coj axe not replicating (senescent). Note that setting 
either of these values to zero results in a uniform rate of birth/ 
death. Effects of varying proliferation/ death parameters are shown 
in Fig. S2. The parameters governing proliferation, in particular b 
and co/,, have much larger influence on the final differentiated cell 
population size than parameters governing death/removal. The 
steepness of the switch does not substantially change the age 
distribution. 



~^ + ~^a =( rJ i<(a)-Lik(a))Pk + i'*Gk-\(a)p k - l . (2c) 

Here u* is the mutation rate per cell per unit time and pi(a,i) is the 
number of progenitor cells of "age" a from the subpopulation with 
i mutations. We assume 0 < i < M, and no back mutation is 
allowed. 

Let fii(t) be the number of stem cells with i mutations at time t. 
Let otDj be the probability of a symmetric division that gives rise to 
two differentiated cells, %a,i be the probability of an asymmetric 
division that gives rise to one stem cell and one differentiated cell, 
and as j be the probability of a symmetric division that gives rise to 
two stem cells. Then 

Pi(0,t) = (2oio,i + KA,i)ni(t). (3) 



Stem Cells 

Constant stem cell population size. To model the evolu- 
tionary dynamics of a stem cell population under strict homeostasis 
(resulting in a fixed stem cell population size), we used the Moran 
stochastic process for M+\ different types, with mutations 
between types and neutral fitness [50]. Let the number of 
individuals carrying each possible number of mutations be given 

by n = (no,ni,...,/iM), where J^ ( ._ 0 =N SC . We considered two 
versions of this model, with and without dedifferentiation. In both 
cases, we assumed that each stem cell divides, on average, every 
T gm chronological time units. Thus, in a population of size N sc , 
the average time between divisions was T gen /N sc . 

In the first model, no dedifferentiation of progenitor cells was 
possible. Every T gm /N sc time units, a single randomly chosen stem 
cell j was removed and one cell / was born with probability given by 



If we neglect mutation, the steady wave-form solutions of Eq. 2 
have the form 



Pi(a,i) = aHj(f — a)e r ^ a \ 



(4) 



where a = 2oq + o^ ,- is the average number of stem cells of type i 
produced per division and r,(a)= — Hi(s))ds is the age- 

dependent growth rate of the differentiated cell population (Text 
SI). Hence, the long-term age distribution is largely determined by 
the functional forms of the differentiated cell birth and death rates 
(Fig. SI and S2.) Altered birth and death rates due to mutations 
can result in mutant subpopulations growing to higher plateaus in 
size, but the final population size will be bounded. Our PDE 
system can be easily modified to have a maximal carrying capacity 
Ki for each sub-population. This does not qualitatively change the 
age distribution of progenitor cells (Fig. SI) and does not 
significandy affect the fraction of (-mutation cells in the total 
progenitor population (Fig. S3), so we do not consider it further. 



P(n^n + e i -e i )=-f 

Nsc 



E 

h = 0 



m,j, 



nh 



(6) 



where my is the probability of changing to type i from type h per 
replication event, and ej is a unit vector with 1 in the y'th column. W e 
considered a linear cascade of mutations in which the muta- 
tion matrix [my] is given by 



(7) 



i-m o o ••• o 

Ml 1 — u 2 0 • • • 0 

0 w 2 0 ••• 0 

0 0 0 ••• 1 



We also considered a version of the model in which dedi- 
fferentiation of two-mutation differentiated cells was allowed, but 
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the total stem cell population size remained fixed. In this model, 
the probability of death of a y'-mutation stem cell and birth of an i- 
mutation cell was given by 



be obtained from the backward Kolmogorov equation for the 
probability of fixation / of a gene with initial frequency po before 
time t: 



P(n->n + ej-ej) = 



(1 



E 



mi,h 



+ <5i,26 



\ a p 2 (a)da 



(8) 



where e is the proportion of cells in the stem cell pool that come 
from dedifferentiated cells at each replication event, and S i2 is the 
Kronecker delta function signifying that that only two-mutation 
progenitor cells can dedifferentiate. Here 



/ = 



l l p 1 (a)da 



is the proportion of two-mutation cells of all ages in the progenitor 
population, given that Pi(a) is the density of differentiated cells of 
age a carrying i mutations. We also considered a version of the 
model in which all progenitor cells, regardless of the number of 
mutations, could dedifferentiate (Text SI). 

Because the Moran model has been studied extensively, we were 
able to use several existing results on the time to emergence and 
fixation of mutants. Let %m be the first time at which an individual 
carrying M mutations emerges who will go on to fix in the 
population. We focus on the case M = 2 because sequencing of acute 
myeloid leukemia genomes suggests that there are 2 driver mutations 
present [54]. (See Discussion for more details.) Using branching 
process approximations, Durrett et al. [34] calculated the waiting 
time for the Moran model under neutral drift of prior mutants. For 
M = 2, the probability density function for t 2 is given by 



\- e - 2 'l x 

<j>(N sc ui ^u 2 p fix x 2 = t) « 1+e _ 2>/i exp 
This simplifies to 



-II 



t x _ e -im 

0 l+ e -2«7A l 



dt . (9) 



(j>(N sc ui y /u 2 p fix T 2 = t)x tanh(?//l)(cosh(//l)) \ (10) 

where k = N sc U\, and pr lx is the probability that a single mutant 
individual will fix in a population of size N sc . For neutral drift, 
Pfix = 1 /Nsc, an d for weak selection 



1 — exp( — 2.?) 
Pfix= l-exp(-2. ? iV sc )' 



(11) 



where M-mutation cells have advantage s« 1 [50]. 

The time to fixation of the subpopulation with M = 2 mutations 
is a sum of two random variables: the time z 2 until appearance of a 
successful two-mutation cell (Eq. (9)) and the waiting time Xf a from 
the time that mutant first appears until that mutant fixes [55]. 
Note that this time is given in units of stem cell generation times 
T gm . The probability density function of the total fixation time, 
Tf lx , is given by the convolution 

t T flx 

<t>(T fix )= (KWTju-Odt (12) 
Jo 

of the probability density functions (j) for time to first appearance of 
successful mutant and <j> for time it takes that mutant to fix. <j> can 



Sf(p 0 J) Po (l-po)d 2 f(p 0 ,t) ,df(p 0 ,t) 
— 2 \- spoil -po) 



8t 



2N SC 



dpo 



(13) 



subject to boundary conditions f(l,t)= 1 and /(0,?) = 0 and initial 
condition /(/?o,0) = S(po). Dividing/ by the ultimate probability of 
fixation Pf lx and differentiating with respect to t, we obtain the 

probability density function for (j> as a function of initial allele 
frequency po [56]. 

Variable stem cell population size. Our previous stem cell 
models couple birth and death events to keep the population size 
fixed, but we next decoupled these events to allow for a 
stochastically varying population size. For clarity, we refer to the 
total stem cell population size in this model by S(t). Again 
assuming that the average replication time of a stem cell is T gen , 
the interval between birth/ death events in this stochastic stem cell 
model corresponds to T gen /S(t) time units in the progenitor cell 
model. We assume that homeostasis in the stem cell pool is 
maintained by control of cell fate upon division, and that each 
stem cell can produce zero, one, or two stem cell offspring. For 
example, the possible offspring from a zero-mutation stem cell are: 
two differentiated zero-mutation cells with probability Op, one 
zero-mutation differentiated cell and one zero-mutation stem cell 
with probability 0^(1— u), one zero-mutation differentiated cell 
and one one-mutation stem cell with probability <1aU, two zero- 
mutation stem cells with probability — 2u), and a one zero- 
mutation stem cell and one one-mutation stem cell with 
probability 2uus. (For simplicity we assume the probability of 
both offspring carrying new mutations to be negligible.) In general, 
a stem cell carrying i mutations can produce k stem cell offspring 
carrying j mutations with probability P'jik) given by 



P\(0) = a B , Pi(l) = a A (l-u) + oc s u, 
P i + 1 (l) = oi A u + asu, P\(2) = a s (l - 2u). 



(14) 



With constant division probabilities old&s an d O-a, this model is a 
rescaled Galton-Watson branching process, and the stem cell 
population either goes extinct in finite time or undergoes expo- 
nential growth when the mean number of stem cell progeny per 
cell, m = 2ols + %a, is greater than one [57] . To describe a stem cell 
population under homeostasis, the probabilities au, %a, and ccs 
must depend on the total stem cell population size S(t). We model 
a carrying capacity K, for stem cells carrying i mutations, such that 

^ =, ' as '' =(1 - ,) xf + i(?'^ =(1 - ,) W- (15) 

As t] approaches one, cca approaches one, so that most of the 
divisions that occur do not change the stem cell population size, 
and we recover the Moran process for all carrying capacity values. 
The parameter £ controls strength of fluctuations about the 
carrying capacities; as £ increases, the fluctuations become smaller. 
In this model, newly emerging mutants can still go extinct due to 
stochasticity, but the total population reaches a quasi-stationary 
regime at population size Kf. Note that, because this is a quasi- 
stationary regime, eventually the stem cell population will go 
extinct, but the expected time to extinction is exponentially 
proportional to Ki for <^ = 1 and 17 = 0 [58]. We chose the carrying 
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capacity to be large enough that extinction of the stem cell 
population does not occur on a physiological timescale, and we 
initialized the stem cell population to be at carrying capacity with 
zero-mutation cells. 

We considered two versions of the variable stem cell population 
size model. In the first case, no dedifferentiation was possible. In 
the second case, differentiated cells with K mutations were allowed 
to dedifferentiate and re-enter the stem cell population. Let S be 
the dedifferentiation rate per two-mutation progenitor cell per unit 
time. Then the mean number of cells dedifferentiating in the 
interval T gen /S(t) between two stem cell replications is 

1=s m m=s W)\ P2(a,t)da - (16) 

To introduce dedifferentiated cells into the stem cell population, at 
each replication event we calculate the mean number of 
dedifferentiated progenitor cells X and update the stem cell 
population: 

[«o,»i ,n 2 ] -» [n 0 ,m ,ri2 + Poissonfl)] , (17) 

After the total stem cell number is updated, the probabilities of 
reproduction are re-calculated using Eq. (14), and reproduction is 
carried out. 

Model Parameters 

Parameters used are summarized in Table 1 . We used parameter 
estimates from the human hematopoietic system because param- 
eters for other cancers are less well known. We used M = 2 as the 
number of necessary mutations to develop a cancerous phenotype. 
Although it has been estimated that for the human hematopoietic 
system there are 1 1,000-22,000 stem cells [59], which give rise to all 
blood and immune system cells, most of these cells are quiescent and 
only divide when body sustains an injury and needs to repopulate 
the hematopoietic system. Our model only considers actively 
dividing stem cells, which have been estimated by various methods 
to number around 100 [32,60]. The entire actively dividing stem 
cell population has previously been modeled as turning over once 
per year [32] , but most recent estimates have an individual stem cell 
dividing every 25—50 weeks [61]. However, this is likely an over- 
estimate, as it is difficult to distinguish between actively dividing and 
quiescent stem cell populations. We assume that an active stem cell 
divides every 20 weeks, which when multiplied by N sc results in an 
active stem cell population turnover time of T gen = 5 weeks. (The 
entire stem cell population including quiescent cells turns over on a 
much longer timescale.) 

Whereas the size of the active hematopoietic stem cell pool is 
small, the number of progenitor cells such as granulocyte, 
erythroid, monocyte, and megakaryocyte colony-forming units 
(CFU-GEMM) and granulocyte and monocyte colony-forming 
units (CFU-GM) is much larger. There are approximately 10 s 
CFU-GEMM cells and 10 8 CFU-GM cells [62]. There are 
estimates that each CFU-GEMM may contribute to hematopoi- 
esis for an average of 60 days (range of 40-340 days) and that it 
replicates at an average rate of once every 50 days (range of 35- 
285 days) [62]. We track the progenitor populations for L = 20 
weeks, and assume that their proliferative potential rapidly drops 
off after 10 weeks. The maximal proliferation and death rates, bi 
and di were chosen so that 100 stem cells results in 10 s — 10 6 
progenitor cells of all ages. 

Not much is known about the selective advantage 5 provided by 
driver mutations for different cancer types, except that it is small 



(r= 1 +j« 1). Unless stated otherwise, we assume neutral fitness in 
the stem cell pool (s = 0) in our stochastic models throughout the 
paper, to focus on the effect of dedifferentiation. We use a range of 
s = 0 ■ ■ ■ 0.4 for the progenitor cells in the deterministic model. 

Mutation estimates per cell division per gene range from about 
10~ 7 in normal cells to 10~ 2 in the case of chromosomal 
instability [63] . (Note that the rate of epigenetic change has been 
estimated to be orders of magnitude higher than that of genetic 
change and could also play a role in cancer initiation [10].) A 
common value used in many mathematical models is a driver 
mutation rate of w=10~ 5 per division, obtained by assuming a 
somatic mutation rate of 10~ 7 per gene, and about 100 genes that 
could be mutated to give same phenotype [40,45]. In normal 
hematopoietic cells the mutation rate has been measured as 
M=10 per division [64]. 

Note that in the stochastic model, which considers every cell 
division, the mutation rate u can be used as is, but using 
chronological time (i.e., weeks or months) means that this value 
should be multiplied by the average number of divisions per unit 
time to obtain w*. (Mutations that speed up the cell cycle will then 
speed up the apparent mutation rate per unit of chronological time 
in our progenitor model.) The expected number of doublings from 
rii stem to />,• progenitor cells is log 2 (/?;/«;) +2 [46], and the total 

number of progenitors cells of type i is />, k caite^ r ^ > ' <la . Using 
values from Table 1, this results in 8 — 10 cell divisions that take 
place over 10 weeks, so u* «w in equations (2). 

Results 

The coupled system of stem cells and progenitor cells undergoing 
mutation and dedifferentiation we modeled is complex. To 
disentangle the effects of different phenomena, we systematically 
built up the model. We first considered the progenitor population 
alone. We then considered the stem cell population alone, in models 
with strict and variable stem cell homeostasis. Finally, we coupled 
the stem and progenitor populations through dedifferentiation. 

Progenitor Population Alone 

We first considered whether mutation and reproduction in the 
progenitor population could by itself generate a sustained 
population of two-mutation cancerous cells. We thus modeled a 
scenario in which no stem cell mutations occur, so the boundary 
condition to the progenitor population system in equations (2) is 
simply (N sc ,0, ■ ■ ■ ,0). Because selection in the progenitor popula- 
tion might favor mutants, we also assumed that progenitor cells 
with / mutations have a proliferation rate &, = (1 + s)'bo (Eq. (5)). 
This yields a steady-state age distribution of normal and mutant 
progenitor cells (Fig. S2). 

Fig. 2 summarizes results for typical parameter values, showing 
that for M = 2 mutant cells to be an appreciable fraction of the 
population, the mutation rate u and proliferative advantage s must 
both be unreasonably high. This is true both if the total progenitor 
population can grow without bound (Fig. 2a) and if its growth is 
restricted (Fig. 2b). Similar findings are obtained if competition 
between progenitor subpopulations is included in the model (Fig. 
S3). Consistent with previous work [26,27,36], these results show 
stem cell dynamics cannot be ignored in considering time to 
carcinogenesis, so we next considered stochastic models of the 
stem cell population. 

Stem Cell Population Alone 

Our first models for stem cells did not incorporate dedifferen- 
tiation, so the dynamics were entirely governed by the stem cells. 
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Table 1. Parameter values used in numerical simulations. 





Parameter 


Meaning 


Value 


L 


maximum progenitor lifespan 


20 weeks [62] 


T gl .„/N sc 


mean time between stem cell replication/death events 


0.05 per week [45] 


N m 


total number of active stem cells at homeostasis 


1 00 [32] 




turnover time of active stem cell population 


0.05 x 100 = 5 weeks (see text) 


K t 


active stem cell population carrying capacity 


N„ [32] 


hi 


maximal proliferation rate of progenitors with i mutations 


1.5+0.2/ per week (see text) 


dt 


maximal death rate of progenitors with i mutations 


1 .0 per week (see text) 


Pb 


steepness of the proliferation switch 


2 (assumed) 




age at which proliferation switches off 


10 weeks (assumed) 


w d 


age at which death switches on 


0 weeks (assumed) 


u 


stem cell mutation rate per replication 


10- 2 [63] to 10- 6 [64] 


it 


effective progenitor mutation rate (per week) 


same as u (see text) 


n 


probability of stem cell asymmetric division 


0 to 1 


e 


replicating stem cell fraction from dedifferentiated progenitors (constant size model) 


0 to 1 


S 


progenitor cell dedifferentiation rate (variable size model) 


0 to 10 per week (assumed) 


lit) 


mean number of dedifferentiated progenitors per stem cell reproduction event (variable size model) 





doi:l 0.1 371 /joumal.pcbi.l 003481 .t001 



In modeling cancer, the time to carcinogenesis can be defined as 
the time for a single M -mutation cell to emerge, the time for M- 
mutation cells to pass some threshold number or fraction, or the 
time for M-mutation cells to fix in the population. If the mutation 
rate is low (such that N sc li<x\), then all three definitions are 
similar, because the time to emergence of a successful M-mutation 
cell is long compared to the time from emergence to fixation. 
However, there is large uncertainty regarding effective mutation 
rates in carcinogenesis (Table 1), so the assumption of low 
mutation rate may not always be valid, and we thus calculated 
times to fixation. 

We began our stem cell modeling by considering fixed 
population size, corresponding to strict homeostasis. In this 
constant N sc case, we could leverage several analytic results, with 
which our simulations agreed well. Fig. 3A shows a typical 
simulation. The full probability density distribution of time to 
fixation is given by Eq. (12) and agrees well with our simulations 
for high mutation rates (Fig. 3B). The time to emergence of a 
successful mutant is of order 1 / (i/N sc u 3 / 2 ) stem-cell generations 
(Eq. (9)). For normal mutation rates ofw=10~ 6 — 10~ 8 per cell 
division, the mean time until emergence of a two-mutation cell is 
10 8 — 10 stem cell generations, which is very long even with a 
short stem cell generation time. 

Because homeostasis is likely imperfect, we also considered a 
stochastically fluctuating stem cell population size. We found that, 
without dedifferentiation, the distributions of times until fixation 
are very similar for models with and without fluctuations in the 
stem cell population size, as long as we condition on non- 
extinction of the stem cell population (Fig. 3B). This is true for a 
wide range of probabilities of asymmetric division r\ and strengths 
of mean reversion £ (Eq. (15)). This agrees with previous findings 
that demographic stochasticity does not alter fixation times of 
neutral mutants in a large population [65], provided that the 
carrying capacities of the mutants are the same. 

Our results suggest that dynamics within either the progenitor 
or stem cell compartments considered separately do not result in 



carcinogenesis in the hematopoietic system on a realistic time- 
scale, provided that cancer-causing mutations occur at normal 
mutation rates, selection advantages relative to wild-type stem cells 
do not appear until M = 2 mutations, and the stem cell population 
size is constant or varies stochastically around a carrying capacity. 
We thus turned our attention to coupled model systems in which 
progenitor cells can dedifferentiate into stem cells. 

Dedifferentiation with Constant Stem Cell Population 
Size 

For the coupled system, we first considered stem cell homeo- 
stasis caused by strict asymmetric division in the stem cell 
population, so the stem cell population size remains fixed. To 
model dedifferentiation in this case, we built off the Moran model 
and assumed that when a stem cell dies and another enters the 
population, the new entrant comes from the two-mutation 
progenitor population with probability equal to £ times the 
proportion of two-mutation cells in the progenitor population. 
Otherwise the new stem cell comes from replication of another 
stem cell. Roughly speaking, in this model the death of a stem cell 
leaves a opening in the niche, which can potentially be filled by a 
dedifferentiated progenitor cell. The number of progenitor cells 
which can successfully dedifferentiate is controlled by the number 
of niche openings (stem cell deaths), not by the absolute number of 
progenitor cells. 

Typical simulation results are shown in Fig. 4A. We found that 
dedifferentiation dramatically shortens the time to fixation of two- 
mutation cells (Fig. 4B). For small dedifferentiation rates £ < 0.05, 
we also saw good agreement between our simulations and a semi- 
analytical approximation for the time to fixation of two-mutation 
cells with selective advantage £ (Eq. (12)). This agreement suggests 
that under strict stem cell homeostasis, dedifferentiation is 
effectively equivalent to a growth advantage for mutant stem cells. 

Distributions of times to fixation of two-mutation stem cells are 
plotted as a function of both dedifferentiation rate £ and mutation 
rate u in Fig. 4C. Dedifferentiation had two major effects in this 
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Figure 2. Steady-state progenitor distributions in the absence of stem cell mutation. (A, B) Fraction of two-mutation cells as a function of 
mutation rate u and proliferative advantage s for (A) unlimited growth, (B) logistic growth for each subpopulation. (C, D): Corresponding total 
population sizes for (C) unlimited growth, (D) logistic growth. Birth/death rates of progenitor cells are given by Eq. (5) with constant death rate ji= 1 
and sigmoidal birth rate with maximal growth rate fto = 2, 6;=(1 for t = 1 ,2. In (B) the carrying capacity used is 

Ni=200N sc ,N2 = 250N sc ,N] = 300N sc .. Other parameters are as in Table 1. For two-mutation cells to reach appreciable levels in this scenario, 
both the mutation rate and the proliferative advantage must be unreasonably large. 
doi:1 0.1 371 /journal.pcbi.1 003481 .g002 



model: increasing the probability that an emergent two-mutation 
stem cell would fix and reducing the time between emergence and 
fixation. Both of these effects act only after a two-mutation cell has 
been generated in the stem cell population. (Recall that, as shown 
in Fig. 2, the mutation rate and selective advantage must be 
unrealistically high for a nontrivial fraction of two-mutation 
progenitor cells to exist in the absence of underlying two-mutation 
stem cells.) For all mutation rates u, the distribution of times to 
fixation was roughly constant for dedifferentiation rates 



g < 1/ N sc , consistent with population genetics theory that 
selection is only effective when the selection coefficient is greater 
than the reciprocal of the effective population size. For small 
mutation rates u, increasing £ beyond this threshold only 
marginally shortened the total time to fixation. This is because 
in this case the total time to fixation is dominated by the time for a 
successful two-mutation cell to emerge, and dedifferentiation only 
reduces this time by a factor of 1 / JPfix (Eq. (9)), where Pf lx is the 
probability of a emergent two-mutation stem cell fixing. Under 




Figure 3. Times to fixation without dedifferentiation. (A) Typical simulation trajectory for constant stem cell population size and mutation rate 
K=0.1. The numbers of zero-, one-, and two-mutation stem cells are shown in blue, red, and green, respectively. The proportion of two-mutation 
cells in the progenitor population is shown in black. (B) Times to fixation for constant and variable stem cell population size models. Histograms of 
waiting times to fixation of two-mutation cells for constant (blue) and variable stem cell population size with high fluctuations (green, c = 1) and low 
fluctuations (red, c = 100). The semi-analytic distribution of waiting times calculated from Eq. (12) is shown in black. In both panels the mutation rate 

K=0.01. 

doi:1 0.1 371 /journal.pcbi.1 003481. g003 
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neutrality pjj x = \/N sc , so for our model with N sc = 100, dediffer- 
entiation can shorten the time to emergence by at most a factor of 
10. The dedifferentiation rate needed to significantly change this 
waiting time scales linearly with N S c (Fig. S4C). Hence, for larger 
stem cell population sizes, a small dedifferentiation rate would have 
a larger effect. For high mutation rates u, the effect of dedifferen- 
tiation is more dramatic, because the time from emergence to 
fixation of two-mutation cells, which dedifferentiation also shortens, 
is comparable to the time to emergence (Fig. 4D). 

The model considered in Fig. 4 assumes that only two-mutation 
progenitor cells can dedifferentiate. We also considered an 
alternate model in which any progenitor cell can dedifferentiate 
(Text SI). In this alternate model, dedifferentiation again had little 
effect for e< \/N sc - P as t that threshold the effect was substantial, 
because in this model dedifferentiation speeds up the time to 
emergence of two-mutation cells, because one-mutation cells fix 
much more quickly when they too can dedifferentiate (Fig. S4D). 
In addition, we considered the case in which the dedifferentiation 
rate is additionally weighted by the progenitor proliferation rate, 
and our results did not change qualitatively (Text SI, Fig. S4B). 

Our analytical and numerical results suggest that, with intact 
homeostasis in the stem cell population and normal mutation 
rates, dedifferentiation plays a fairly minor role in speeding up the 
time to cancer initiation. We thus turned to consider the case in 
which homeostasis is not strict. 

Dedifferentiation with Variable Stem Cell Population Size 

In the previous section, we assumed that the stem cell 
population size was constant because homeostasis was maintained 



by all divisions being strictly asymmetric. Consequently, dedif- 
ferentiated progenitor cells could only occupy newly created 
openings in the stem-cell niche created by a death event in the 
stem cell population. Because homeostasis is likely maintained at 
the population level [66], with each stem cell division producing 
not stricdy one stem cell but rather on average one stem cell, we 
next considered a model in which the stem cell population could 
stochastically fluctuate around a carrying capacity. In this model, 
stem cell homeostasis was maintained by dynamically altering the 
probabilities of the three possible outcomes of a stem cell division: 
two stem cells, one stem and one progenitor cell, or two 
progenitor cells (Eq. (15)). Two-mutation progenitor cells each 
had a probability per unit time of dedifferentiating, and 
dedifferentiated cells were simply added to the stem cell pool. 
Thus in this model the total influx of dedifferentiated cells 
depended on the total number of two-mutation progenitor cells, 
not on the creation of openings in the stem cell niche. (Note that, 
in our previous model with constant stem cell population size, the 
rate of dedifferentiation per reproduction event was denoted e. 
To distinguish the present model, we denoted the progenitor 
dedifferentiation rate per cell per unit time as 5.) Again, we asked 
whether dedifferentiation substantially speeds the time to 
carcinogenesis. 

Fig. 5A and 5B show typical results from this model for a 
moderate dedifferentiation rate S. After a waiting time, the 
population of stems cells began to grow exponentially, because the 
influx of dedifferentiated two-mutation progenitor cells exceeded 
the capacity of stem-cell division homeostasis. For larger 
dedifferentiation rates, the exponential growth rate is larger 




Figure 4. Times to fixation with dedifferentiation for constant stem cell population size. (A) Typical simulation trajectory with 
dedifferentiation (e = 0.02) for the same random number seed as Fig. 3A. Blue: zero-mutation stem cells, Red: one-mutation stem cells, Green: two- 
mutation stem cells. Black: proportion of two-mutation cells in the progenitor population. (B) Distributions of times to fixation of two-mutation cells 
under strict stem cell homeostasis. Normalized histograms (dots) and analytical approximations (solid lines) are shown for » = 0.01 for zero 
dedifferentiation (red; c = 0) and non-zero dedifferentiation (black; c = 0.02). (C) Median times to fixation of two-mutation cells (solid lines) and inter- 
quantile ranges (shaded regions) versus dedifferentiation rate c and mutation rate u. (D) Mean times to emergence of a successful two-mutation stem 
cell (solid lines, Eq. (9)) and fixation of such cells (dotted lines, Eq. (12)) in Moran models with selection coefficient s. Black curve indicates first 
appearance of two-mutant cell. 
doi:1 0.1 371 /journal.pcbi.1 003481. g004 
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(Fig. 5C and 5D), and the distribution of progenitor ages can be 
distorted, with many young cells, as seen in Fig. 5E and 5F. 

Exponential growth eventually occurs whenever the dediffer- 
entiation rate exceeds a threshold <5 OT( . Solving self-consistently for 
the influx of dedifferentiated cells and the growth rates of the stem 
and progenitor cell populations, we obtained an integral equation 
for the growth k 

k = asf e- ka <? (a) da- X -^p-, (18) 
J o T gen 

which provides an excellent fit to the numerical simulations (Fig. 5 
and 6A,B). (For derivation details, see Text SI.) Setting this growth 
rate k to zero, we found 

a " aT gen J 0 °° e'^da 

Here rj is probability of asymmetric stem cell division (producing 
one stem and one progenitor cell), and T gm is the mean time 
between stem cell divisions. (Note that if r\ = 1, this model reduces 
to the Moran model with the population size monotonically 
increasing due to dedifferentiation.) Lasdy, a = 20(^2 + 04,2 in Eq. 
(19) is the average number of progenitor offspring produced by a 
two-mutation stem cell. Because ocd^ changes as the system 
attempts to maintain stem-cell homeostasis, a is actually a 
stochastic variable that depends on the stem cell population size. 
During exponential growth axl — rj, because the probability of 
symmetric divisions that give rise to two stem cells goes to zero, 
and all new stem cell growth comes from dedifferentiated 
progenitor cells. In Eq. (19), r(a) is the growth rate of two- 
mutation progenitor cells as a function of age a, so Jg™ e'^da is the 
number of progenitors produced by one two-mutation stem cell. 
Increasing the amplification of mutant stem cells into progenitors 
increases the net dedifferentiation rate, lowering the threshold 
5 crit . Because the threshold 5 crit depends on the age distribution of 
the two-mutation cells, for a given (small) rate of dedifferentiation 
<5, evolving a mutant that proliferates faster (increasing e^"') can 
destabilize a system in which the number of cancerous cells is 
stable and take it into exponential growth regime. 

The dependence of the critical dedifferentiation rate 8 crU on the 
growth-rate advantage s of two-mutation progenitor cells and 
probability rj of asymmetric cell division is shown in Fig. 6B. The 
critical <5 decreases rapidly as the selective advantage of two- 
mutation cells increases. Increasing rj or T ge „ also lowers the 
critical dedifferentiation rate, because homeostasis is less effective 
when asymmetric stem cell divisions are less frequent. Note that 
the exponential growth rate k does not depend on the mutation 
rate (Fig. S5 A), and although the critical S given by ( 1 9) needed for 
exponential growth is a function of the probability of asymmetric 
division rj, the actual growth rate k and the time to exponential 
growth are not significantly affected by changing rj (see Fig. S5B). 

For dedifferentiation rates <5 below S crit , two-mutation stem cells 
eventually fix in the population, but for S>S crit , the stem cell 
population is likely to begin exponential growth before fixation of 
two-mutation stem cells. Thus in Fig. 6C and 6D we report the 
time to carcinogenesis as the time for the two-mutation stem cell 
population to exceed N sc , the nominal carrying capacity of the 
stem cell compartment. In this case of stochastic stem cell 
homeostasis, dedifferentiation can dramatically shorten the time to 
carcinogenesis, even for low mutation rates u. This is because the 
first two-mutation stem cell often arises not from direct mutation 
of a stem cell, but rather from dedifferentiation of a progenitor cell 



generated by mutations within the progenitor compartment 
(Fig. 6E). Although mutations in the progenitor compartment do 
not affect a large fraction of progenitors, because the number of 
progenitor cells is so large, the absolute number of two-mutation 
progenitor cells is non-negligible. Thus even small rates of 
dedifferentiation can have dramatic effects. This is in contrast to 
the case of strict stem cell homeostasis, in which the absolute 
number of two-mutation progenitor cells was unimportant, 
because they needed an opening in the stem cell niche to 
successfully dedifferentiate. 

Our results show that the case of stochastically controlled stem 
cell homeostasis is qualitatively different from the case of strict 
homeostasis. If homeostasis is controlled at the population level 
(where stem cell decisions between symmetric and asymmetric 
division are stochastic), dedifferentiation can overwhelm it, leading 
to exponential growth of the stem cell population. Moreover, if 
dedifferentiated cells do not depend on openings to colonize the 
stem cell niche, dedifferentiation can dramatically hasten the time 
to carcinogenesis, even for low mutation rates. 

Discussion 

Progression to cancer is associated with expansion of the cancer 
stem cell (CSC) population, but the origin of these CSCs remains 
unclear. Although CSCs may arise directly from adult stem cells, 
they may also arise from somewhat differentiated cells that have 
dedifferentiated and acquired stem cell-like characteristics 
[13,14,18,19,67]. Stems cells replicate indefinitely, giving them a 
long time to accumulate the mutations that drive carcinogenesis, 
but the population of actively dividing stems cells (N sc ) is small. 
Progenitor cells replicate only a small number of times, but the 
population of progenitor cells is typically several orders of 
magnitude larger than the stem cell population. Thus, as a 
population, progenitors undergo many more divisions, potentially 
letting some of these cells acquire mutations that enable them to 
dedifferentiate and drive carcinogenesis. Here, using mathematical 
modeling, we have shown that even a small rate of dedifferenti- 
ation may drastically shorten the time to cancer emergence, even 
for low mutation rates. 

Recent studies suggest stem cell dynamics during homeostasis 
are governed by neutral competition and genetic drift [10,29,30]. 
Traditionally, stem cells were thought to always undergo 
asymmetric division, always yielding a stem cell and a progenitor 
cell, resulting in a fixed stem cell population size. This scenario is 
represented by our first model for stem cell dynamics, based on the 
popular Moran model. It has been recently shown, however, that 
symmetric divisions also occur in adult stem cells and may be the 
predominant form of division [68,69]. Moreover, cancer stem cells 
have been shown to undergo more symmetric divisions than 
normal stem cells [70]. Little is known, however, about how the 
stem cell population size is regulated [29]. Hence, in our second 
model for stem cell dynamics, we made the simplifying assumption 
of an a priori carrying capacity Kj. We considered a density- 
dependent stochastic process, in which the degree of mean 
reversion is controlled through the probabilities of producing zero, 
one, or two stem cell offspring. In this model, the non-constant 
stem cell population size S(t) tends to return to the carrying 
capacity Kj, because the mean number of stem cells produced per 
division is greater than one when S(t) < Ki and less than one when 
S(t)>K,. (Although the stem cell population size could, in 
principle, be maintained by regulating apoptosis rather than 
biasing division, previous modeling suggests that regulating 
division probabilities rather than cell cycle time or removal is 
more important for maintaining homeostasis [46,71].) 
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Figure 5. Exponential growth given varying stem cell population size and dedifferentiation. Total number of stem cells (A, C, E) and 
corresponding final progenitor age distributions (B, D, F) are shown in black. Also shown are best exponential fits of the growth rate (blue) and our 
semi-analytic solution given by Eq. (18) (green). In all panels the probability of asymmetric stem cell division 17 = 0 and the mean reversion parameter 
is {=1. 

doi:1 0.1 371 /journal.pcbi.1 003481 .g005 



If the stem cell population size varies and is regulated by biasing 
division, we found two distinct regimes. If the dedifferentiation 
rate is much less than a critical value, then the initial two-mutation 
stem cell often arises from a normal stem cell, so the time to 
fixation of such a cell is similar to the case with constant 
population size. If the dedifferentiation rate exceeds the critical 
value, however, then the initial two-mutation stem cell often arises 
from a dedifferentiated progenitor cell, so the time to fixation is 
dramatically shorter than the case with constant population size. 
Moreover, in this regime the stem cell population eventually grows 
exponentially, as dedifferentiating progenitor cells overwhelm stem 
cell homeostasis. Note that the threshold between these two 
regimes is independent of the overall mutation rate, if stem and 
progenitor cell mutation rates are proportional. 

When the stem cell population size is constant, dedifferentiation 
simply acts like a selective advantage for mutant stem cells. When 
the stem cell population size is allowed to vary, however, 
dedifferentiation can additionally drive exponential growth of 
the stem cell population. If the stem cell population size N sc is 



constant, our results imply that stem cell dynamics in the coupled 
stem cell-progenitor system can be approximated by a population 
genetics model of the stem cells alone, as long as that model 
includes positive selection. In this case, we found that the 
dedifferentiation rate £ must exceed 1 /N sc to substantially shorten 
the time to cancer acquisition, similar to classical population 
genetics results that the selection coefficient must exceed the 
inverse population size to be effective. Hence, our model predicts 
that in tissues where the niche contains fewer cells, smaller rates of 
dedifferentiation are sufficient to influence the time to cancer. It is 
interesting to note that cancers where dedifferentiation has been 
shown to occur have a small niche size (i.e. intestinal crypts) 
[15,16]. For the hematopoietic system, based on the available 
literature, we assumed that the number of actively dividing stem 
cells is N sc = 100, so the dedifferentiation rate must be «0.01 or 
higher to significantly shorten the time to cancer. 

Here we focus on the hematopoietic system, in which the stem 
cell compartment consists of N sc xl00 active cells, and two 
mutations are necessary for carcinogenesis. For some other 
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Figure 6. Fixation and exponential growth of two-mutation cells with dedifferentiation for variable stem cell population size. A: 

Observed growth rate k of the stem cell population (black curve) and the semi-analytic approximation Eq. (18) (green) for 17 = 0, c = 1, and m = 0.01. 
The vertical line denotes <5„„. B: Analytically predicted critical dedifferentiation rate Sou as a function of asymmetric division probability 17 and the 
growth advantage s of the two-mutation progenitor population. Exponential growth occurs for <5><5 tT , ( . C: Normalized histogram (stars) of waiting 
times for exponential growth of the stem cell population with stochastic homeostasis and dedifferentiation for « = 0.01, 5 = 0.01. For comparison the 
histogram (red and black dots) as well as the analytical distributions of times to fixation given strict homeostasis for £ = 0.01 and e = 0 are also shown. 
D: The median and inter-quantile range of times to first occurrence of N sc = 100 two-mutation stem cells, given stochastic homeostasis and a range of 
dedifferentiation rates S. For comparison, the waiting times to fixation for N se = 100 given strict homeostasis (shaded areas) for the equivalent value 
of e are also shown. E: The probability that the first two-mutation stem cell arose from mutation in the stem cell compartment, rather than 
dedifferentiation. Vertical line denotes S cri ,. Parameters for all simulations given in Table 1. 
doi:1 0.1 371 /journal.pcbi.1 003481. g006 



cancers, such as colon cancer, the number of stem cells per 
compartment is much smaller, there are many compartments, and 
the number of necessary mutations is larger. For high mutation 
rate, the mean time to fixation scales linearly with N sc (see Fig. 
S4C). So in cancers with small N sc two-mutation stem cells will fix 
much faster. However, the need to accumulate more mutations 
will slow carcinogenesis. We expect, however, that the qualitative 
effects of dedifferentiation will be similar to the hematopoietic 
system we analyzed. 

The fact that mutants take a long time to reach an appreciable 
fraction of the stem cell population is not typically considered in 
the cancer modeling literature, which often makes an implicit 
assumption that a newly emerged mutant cell will not go extinct 
and will fix quickly. Our results show that, for high mutation rate, 



the time for a mutation to fix in the population is comparable to 
time for a successful mutant to first emerge, in accordance with 
classical results of Kimura and Ohta [56]. This is especially 
important if division events are rare and the population size is 
large. Considering the time to some predetermined diagnosis 
threshold is similar to considering the time to fixation, because the 
time between a selected mutation becoming common and fixing is 
typically short [72]. Hence, elevated mutation rate (genetic 
instability) may not speed up time to carcinogenesis as much as 
is typically assumed, suggesting that some form of selection 
(potentially through dedifferentiation) is necessary. Most tumors 
accumulate hundreds of mutations, but the number of necessary 
"driver" mutations depends on the type of cancer. We considered 
M = 2 mutations, because sequencing of acute myeloid leukemia 
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genomes suggest that there are two driver mutations present [54] . 
Moreover, recent findings on induced pluripotent stem cells also 
suggest M = 2, as loss of both copies of the tumor suppressor 
protein p53 [73] or the activation of two oncogenes [67] may be 
necessary for dedifferentiation. Disabling both copies of p53 
improves the efficiency of reprogramming to a stem-like state and 
greatly enhances the production of induced pluripotent stem cells 
[74,75]. The loss of p53 also leads to the emergence of tumor cells 
bearing functional and molecular similarities to stem cells [22,73]. 
Finally, inactivation of p53 changes the ratio of symmetric to 
asymmetric division in mammary stem cells, allowing the total 
stem cell population to escape homeostasis [70]. 

Our model only considers actively dividing stem cells, which in 
the human hematopoietic system have been estimated to be 
roughly 100 [32] out of 11,000-22,000 total stem cells [59]. A 
more complete model would consider both the active and 
quiescent stem cell populations. Transitions between these states 
may be influenced by the progenitor population size, potentially 
acting as a negative feedback and regulating the proliferation of 
cancer stem cells. In our models, cancerous cells take over the stem 
cell population, but the ratio of cancer progenitor cells to cancer 
stem cells is fixed by the progenitor growth process. Even when 
dedifferentiation drives exponential growth of the stem cells, it is 
their absolute number that increases, not their proportion in the 
population. This is in concordance with some in vitro studies, which 
suggest a fixed proportion of CSCs in a tumor [18]. 

Many theoretical models find that in order to accumulate 
multiple mutations on a reasonable time scale, the onset of 
elevated mutation rate (i.e., genetic instability) should be an early 
event in tumorigenesis (reviewed in [4,7]). The importance of 
genetic instability, however, depends on assumptions about 
symmetric self-renewal and differentiation of stem and progenitor 
cells. In particular, mutations that alter stem cell division or make 
committed progenitors somewhat immortal may also lead to an 
early onset of cancer, diminishing the impact of genetic instability 
[45]. Similarly, our results show that different assumptions about 
how dedifferentiation occurs (frequency-dependent reproduction 
versus absolute numbers of dedifferentiating cells) dramatically 
alter time to carcinogenesis. 

A large body of modeling work in this area (reviewed in section 
Prior Mathematical Models) has focused on calculating the time to 
carcinogenesis under the assumption of constant population size 
(not specifying the mechanism of homeostatic regulation). We 
compared the times to multiple mutation acquisition in our 
constant and variable stem cell population size models and found 
that without dedifferentiation both models yield similar results. 
With dedifferentiation, however, we found that the two models 
differ substantially. We explicitly considered different ways that 
homeostasis can be maintained in the stem cell population, and 
showed that these assumptions can lead to very different results. 
Our results suggest that if homeostasis is controlled through 
division asymmetry and if de-differentiated cells do not depend on 
openings to colonize the stem cell niche, then for de-differentiation 
rate larger than a critical threshold, the cancer stem cell will most 
likely originate in a progenitor cell that has undergone de- 
differentiation. This is a prediction of our model that can be 
experimentally tested using inducible genetic labeling, the same 
technique that permitted lineage-tracing experiments allowing 
quantification of symmetric versus asymmetric divisions [29,30]. A 
similar method was previously used to identify oligodendrocyte 
precursor cells as the tumor cell of origin in glioma [76] . 

Our model contributes to the existing literature on the trade-offs 
between symmetric and asymmetric divisions of the stem cell 
population in stem cell-driven [35,51,77]. To our knowledge, our 



model is the first to quantify the effects of dedifferentiation on the 
time to carcinogenesis. There are a number of important aspects 
of homeostasis our model does not consider, such as spatial aspects 
of stem cell position in the niche (we assume the cells are well- 
mixed) or negative feedback to stem cell divisions from the 
progenitor population (we assume progenitors only influence the 
stem cells by dedifferentiating). Lander et al. have previously 
shown that negative feedback control is needed for homeostasis, 
and that feedback regulating replication probabilities is more 
effective than feedback regulating cell cycle lengths [71]. In 
ongoing work, we are investigating what effect including spatial 
structure and feedback from progeny will have on dedifferentiation 
times. The effect of spatial structure on mutation acquisition is still 
not fully resolved. Some groups argue that time to acquire M = 2 
mutations is actually decreased in a spatial model compared to the 
space-free model [78]. Other groups find that time to multiple 
mutations is increased when space is considered [79,80]. Like 
other mathematical models, our model suggests that eradication of 
cancer is dependent on eradication of cancer stem cells [81-83]. 
The potential for progenitor cells to dedifferentiate and repopulate 
the stem cell compartment, however, may complicate successful 
treatment. Our work suggests that further progress in understand- 
ing initiation and treatment of cancer requires a more detailed 
understanding dedifferentiation and of stem cell homeostasis. 

Supporting Information 

Figure SI Progenitor cell dynamics for a constant stem 
cell population. Fixed number of zero-mutation stem cells 
A^, f = 100. (A) Proliferation and death rates ff(a) (red curve) and 
[i(a) (blue curve) used in the simulation and the resulting age- 
dependent growth rate r(a) (black curve). (B) The total number of 
progenitor cells given by Eq. (1) (solid line) and Eq. (S6) (dotted 
line) as a function of time. (C) The steady state age distribution 
with no competition between progenitor cells is given by Eq. (1) 
(solid line) and with logistic growth by Eq. (S6) (dotted line). (D-F) 
Progenitor dynamics for for K = 3 subpopulations with mutation 
rate u = 0.001 . (D) Age-dependent birth rates a(a) given by Eq. (5). 
The death rate is constant. (E) The total number of progenitor cells 
as a function of time for zero-mutation (blue), one-mutation (red), 
and two-mutation (green) subpopulations. (F) The steady state age 
distribution for zero-mutation (blue), one-mutation (red), and two- 
mutation (green) progenitor cells. 
(EPS) 

Figure S2 Robustness to parameter variations in pro- 
liferation/death rates in progenitor model. The effect of 
parameters in age-dependent birth and death rates, a(a) and fi(a) 
(given by Eq. (5)) on the age-dependent growth rate 
r(a) = Jq (ct(.s) — fi(sj)ds and the steady state age distribution e r ' a ' 
in our models. (A,B) Effect of varying maximal growth rate b 
between 0 and 2. (C,D). Effect of varying maximal death rate d 
between 0 and 5. (E,F) Effect of varying the location (age of onset) 
of the proliferation switch co/, between 0 and 5. (G,H) Effect of 
varying the age at which the apoptosis switch a>d is turned on 
between 5 and 10. (I J) Effect of varying the steepness of the 
proliferation switch p b between 0 and 5. (K,L) Effect of varying the 
steepness of the apoptosis switch p d between 0 and 5. 
(EPS) 

Figure S3 Steady-state progenitor distributions in the 
absence of stem cell mutation but with progenitor 
competition. Top: The fraction of mutant cells as a function 
of mutation rate u and proliferative advantage s for (A,C) local 
(age-dependent) competition between subpopulations given by Eq. 
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(S7), and (B,D) global competition between subpopulations given 
by Eq. (S8). Bottom: Corresponding plots of total cell density. 
Basal dynamics are constant death rate fi = 1 and sigmoidal birth 
rate with maximal growth rate bo = 2, &; = (1 +s)bi-i for i= 1,2. 
The same carrying capacity is used for all simulations: 
N\ = 200N SC , N 2 = 250N SC , N 3 = 300N SC . Note that there is a 
sharp transition zone at which mutant cells go from nearly zero 
fraction of total population to majority of the differentiating cell 
population. However, the mutation rate u and proliferative 
advantage s at which this is observed is unreasonably high, just 
as for the model without progenitor competition (Fig. 2). 
(TIF) 

Figure S4 Comparison of two Model I variants with all- 
mutant progenitor dedifferentiation and two-mutant 
progenitor dedifferentiation. (A) Fixation time distributions 
in constant stem cell population size model for potential 
dedifferentiation of only two-mutation progenitors (red, Eq. (8)) 
and potential dedifferentiation of all progenitor cells (blue, Eq. 
(S9)). (B) Fixation time distributions in constant stem cell 
population size model with dedifferentiation of all progenitor 
cells. Blue: all progenitor cells equally likely to dedifferentiate with 
dedifferentiation probabilities given by Eq. (S9). Red: all 
progenitor cells can dedifferentiate with dedifferentiation proba- 
bility weighed by birth rate given by Eq. (S10). Progenitor 
dynamics without competition (Eq. (2)). Green: all progenitor cells 
can dedifferentiate with dedifferentiation probability weighed by 
birth rate given by Eq. (S10). Progenitor dynamics with local 
competition given by Eq. (S7). Dedifferentiation rate used is 
s = 0.02, mutation rate is w = 0.01. (C) Mean + standard deviation 
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